Automatic Detection Of Fires On Earth&#39;s Surface And Of Atmospheric Phenomena Such As Clouds, Veils, Fog Or The Like, Using A Satellite System

ABSTRACT

A method for automatically detecting fires on Earth&#39;s surface using a satellite system is provided. The method includes acquiring multi-spectral images of the Earth at different times, using a multi-spectral satellite sensor, each multi-spectral image being a collection of single-spectral images each associated with a respective wavelength (λ), and each single-spectral image being made up of pixels each indicative of a spectral radiance (R λ ) from a respective area of the Earth. The method also includes providing a model relating the spectral radiances (R λ ) of each pixel in multi-spectral images acquired at different times and physical quantities representing thermodynamic phenomena occurring on the Earth&#39;s surface, in the Earth&#39;s atmosphere, and related to the Earth and the Sun relative positions. The method further includes computing for each pixel, at a given time, at least the physical quantity in the model representing a possible fire on the Earth&#39;s surface.

TECHNICAL FIELD OF THE INVENTION

The present invention concerns with automatic detection of fires on Earth's surface and of atmospheric phenomena such as clouds, veils, fog or the like, by means of a satellite system, in particular by exploiting multi-spectral data acquired by multi-spectral sensors of a geostationary or polar satellite system.

BACKGROUND ART

As is known, multi-spectral images are images acquired by Remote Sensing (RS) radiometers, each acquiring a digital image (in remote sensing, called a scene) in a small band of visible spectra, ranging from 0.4 μm to 0.7 μm, called red-green-blue (RGB) region, and going to infra-red wavelengths of 0.7 μm to 10 or more μm, classified as NIR (Near InfraRed), MIR (Middle InfraRed), FIR (Far InfraRed) or TIR (Thermal InfraRed). A multi-spectral image is hence a collection of several single-spectral (single-band or monochrome) images of the same scene, each taken with a sensor sensitive to a different wavelength.

Different fire detection techniques, based on threshold criteria and contextual algorithms, have been developed for multi-spectral polar sensors and, in the last years, for geostationary sensors. For a detailed discussion of these techniques reference may, for example, be made to Kaufman, Y. J., Justice, C. O., Flynn, L. P. Kendal, J. D., Prins, E. M., Giglio, L. Ward, D. E. Menzel, W. P. and Setzer, A. W., 1998, Potential global fire monitoring from EOS-MODIS, Journal of Geophysical Research, 103, 32215-32238, and Giglio, L., Descloitres, J., Justice, C. O. & Kaufman, Y. J. (2003), An enhanced contextual fire detection algorithm for MODIS, Rem. Sen. Environment, 87:273-282.

Multi-spectral sensors on polar satellites are characterized by a relatively high spatial resolution, but, due to the long revisit time of polar satellites, the promptness needed for effective fire detection purposes cannot be achieved, even combining all existing multi-spectral polar sensors. On the contrary, multi-spectral geostationary sensors provide very frequent acquisitions, e.g. every 15 minutes for the MSG SEVIRI (Spinning Enhanced Visible and Infra Red Imager) sensor, though characterized by a lower spatial resolution (3×3 km² and above for infrared channels), which could prevent small fires from being detected.

In order to overcome the spatial resolution limitations, a physical model-based approach for sub-pixel fires detection from geostationary sensors data was recently proposed by E. Cisbani, A. Bartoloni, M. Marchese, G. Elisei, A. Salvati, Early fire detection system based on multi-temporal images of geostationary and polar satellites, IGARSS 2002, Toronto, 2002, and Calle, A., Casanova, J. L., Moclan, C., Romo, A. J., Costantini, M., Cisbani, E., Zavagli, M., Greco, B., Latest Algorithms and Scientific Developments for Forest Fire Detection and Monitoring Using MSG/SEVIRI and MODIS Sensors, IEEE, 2005, 118-123.

In particular, an analytic Radiative Transfer Model (RTM) is proposed which characterizes the radiative phenomena that determine the sensor-detected energy, expressed by means of radiances R_(λ) (W/m²/sr/μm) for each band λ in atmospherically transparent windows in Near Infrared (NIR), Middle Infrared (MIR) and Thermal Infrared (TIR) spectral regions. As shown in FIG. 1( a), the radiance R_(λ) collected by a remote satellite sensor is the sum of the solar radiance R_(S,λ) reflected by the ground, the atmospheric thermal radiance R_(A,λ) (both the up-welling and the down-welling components), and, finally, the thermal emission of the ground. Given the background temperature T_(B), the emissivity ε_(λ) of the Earth's surface, and the transmittance of the atmosphere τ_(λ) between the Earth's surface and the sensor, the RTM can be expressed as:

R _(λ)=ε_(λ)τ_(λ) B _(λ)(T _(B))+R _(A;λ) +R _(S;λ),  (1)

where, B_(λ)(T) is the Planck black-body emission at temperature T and wavelength λ. Other RTM models can be exploited as well.

According to C. C. Borel, W. B. Clodius, J. J. Szymanski and J. P. Theiler, Comparing Robust and Physics-Based Sea Surface Temperature Retrievals for High Resolution, Multi-Spectral Thermal Sensors Using one or Multiple Looks, Proc. of the SPIE '99, Conf. 3717-09, the main contribution to the transmittance τ_(λ) in the atmospheric windows in the NIR and TIR regions comes from the atmospheric water vapor content and the relations between transmittance and water vapor can be quite appropriately be parameterized by the following expression:

$\begin{matrix} {{\tau_{\lambda}\left( {W,\vartheta} \right)} = {\exp \left( {- \left\lbrack {\frac{A_{\lambda}}{\cos \; \vartheta} + {B_{\lambda}\left( \frac{W}{\cos \; \vartheta} \right)}^{C_{\lambda}}} \right\rbrack} \right)}} & (2) \end{matrix}$

where W is the total water vapour along the path ending/starting at/from the examined pixel and having a zenith angle θ. Parameters A_(λ), B_(λ) and C_(λ) depend (at least) on the wavelength λ and can be estimated via several MODTRAN (MODerate resolution atmospheric TRANsmission) simulations (computer program designed to model atmospheric propagation of electromagnetic radiation from 100-50000 cm⁻¹ with a spectral resolution of 1 cm⁻¹) and regression methods. Other models/methods to estimate τ_(λ) can be considered.

The water vapour W content can be estimated as described in Eumetsat Satellite Application Facility, Software Users Manual of the SAFNWC/MSG: Scientific part for the PGE06, SAF/NWC/INM/SCI/SUM/06, issue 1.0, January 2002, but other methods can be considered.

The solar term R_(S,λ) can be calculated as described in the aforementioned Potential global fire monitoring from EOS-MODIS:

$\begin{matrix} {R_{S,\lambda} = {E_{S,\lambda}{\cos \left( \vartheta_{S} \right)}{\tau_{\lambda}\left( z_{SE} \right)}\frac{1 - ɛ_{\lambda}}{\pi}{\tau_{\lambda}\left( z_{ED} \right)}}} & (3) \end{matrix}$

where the E_(S,λ) is the Sun radiance at the top of the atmosphere, τ_(λ) (z_(SE)) is the transmittance along the path between Sun and Earth's surface, τ_(λ) (z_(ED)) is the transmittance along the path between Earth's surface and satellite sensor, and ε_(λ) is the emissivity of the Earth's surface. Other models/methods can be exploited to calculate R_(S,λ).

The atmospheric radiance contribution R_(A,λ) describes a complex phenomenon, characterized by smoke, aerosol, and local atmospheric temperatures hard to be modelled. A possible model is the following:

$\begin{matrix} {R_{A,\lambda} = {\underset{{down} - {welling}}{\underset{}{\left( {1 - ɛ_{\lambda}} \right){{\tau_{\lambda}\left( z_{ED} \right)}\left\lbrack {1 - {\tau_{\lambda}\left( z_{TOA} \right)}} \right\rbrack}{B_{\lambda}\left( T_{A} \right)}}} + \underset{{up} - {welling}}{\underset{}{\left\lbrack {1 - {\tau_{\lambda}\left( z_{ED} \right)}} \right\rbrack {B_{\lambda}\left( T_{A} \right)}}}}} & (4) \end{matrix}$

where, referring to FIG. 1( b), the τ_(λ) (z_(TOA)) is the transmittance along the vertical path between the Earth's surface and the top of the atmosphere, τ_(λ) (z_(ED)) is the transmittance along the path between the Earth's surface and the satellite sensor, and ε_(λ), as in the foregoing, is the ground emissivity.

The Dozier formulation as described in J. Dozier, A Method for satellite identification of surface temperature fields of subpixel resolution, Remote Sensing of Environment, 11 (1981) 221-229 and applied to equation (1) (or to other RTM models) makes a sub-pixel description of the radiative process possible, considering the fire extension (pixel fraction f of a pixel of radiance acquired by the satellite sensor) and fire temperature T_(F):

R _(λ)=ε_(F;λ)·τ_(λ) ·B _(λ)(T _(F))·f+ε _(B;λ)·τ_(λ) ·B _(λ)(T _(B))·(1−f)+R _(A;λ) +R _(S;λ),  (5)

where ε_(F,λ) and ε_(B,λ) are the fire and the background emissivities, respectively, at the wavelength λ.

According to the aforementioned Early fire detection system based on multi-temporal images of geostationary and polar satellites, if two successive acquisitions are considered, the Dozier formulation (5) can be written as follows:

ΔR_(λ,y) ≡R _(λ,t) −R _(λ,t-Δt)=ε_(λ)·τ_(λ,t) ·[B _(λ)(T _(F))−B _(λ)(T _(B))]·Δf  (6)

where t and t−Δt denote two close acquisition times, and Δf=f_(t)−f_(t-Δt), and where the following assumptions are done:

-   -   the background temperature T_(B) is constant between two         consecutive acquisitions (within 15 minutes for MSG SEVIRI         data);     -   the fire emissivity is the same as the non-burning surface         emissivity; and

the solar and the atmospheric terms (R_(S,λ), R_(A,λ)) are considered constant between two consecutive acquisitions (within 15 minutes for MSG SEVIRI data).

The introduced RTM (1) and equations (2), (3), (4), (5), (6) are reliable only if no clouds are in the analyzed scene. Therefore, a reliable cloud masking procedure is necessary to identify the acquisitions that are compatible with the physical model assumptions. Many techniques have been developed for cloud masking by exploiting polar or geostationary sensors. Basically, all these techniques are based on the application of threshold criteria to analytic relations among the different bands of a single acquisition. Methods to retrieve such relations among the spectral bands can be based on physical models or on learning methods based on neural networks, Bayesian networks, support vector machines, all of which require a pre-processing phase for the system training. Also contextual techniques that exploit the spatial information are known in literature.

OBJECT AND SUMMARY OF THE INVENTION

The Applicant has noticed that the analysis disclosed in the aforementioned Early fire detection system based on multi-temporal images of geostationary and polar satellites and represented by the equation (6), where two successive acquisitions are exploited, suffers from rough approximations and strong restrictions which lead to unsatisfactory results for reliable applications. In particular, while the atmospheric contribution in the RTM can be considered slowly changing with time, the estimation of the atmospheric contribution R_(A,λ) given by equation (4) in combination with equation (2) in the RTM is affected by a large uncertainty due to noise and model inaccuracy, and the non-linear form of the model given by equation (5) amplifies this uncertainty, so making the estimation of the fire parameters unreliable.

Moreover, the Applicant has also noticed that, at each new acquisition, the equation (6) is solved and Δf, T_(F) and T_(B) are estimated without exploiting previously computed results, i.e. Δf, T_(F) and T_(B) computed for the previous acquisitions.

The objective of present invention is to provide an improved automatic technique for fire detection from geostationary satellite multi-spectral data which overcomes the limitations of the known techniques.

This objective is achieved by the present invention in that it relates to a method for automatically detecting fires on Earth's surface and atmospheric phenomena such as clouds, veils, fog or the like, by means of a satellite system, as defined in the appended claims.

The present invention stems from the observation by the Applicant that many physical quantities in equation (1), such as the atmospheric contribution R_(A,λ), the solar contribution R_(S,λ), the transmittances τ_(λ) and the emissivities ε_(λ) are highly temporally correlated. For this reason, knowledge and exploitation of many acquisitions of the same scene allows information of interest, such as fires, clouds and physical parameters estimations, to be retrieved with an higher accuracy and robustness than considering only one or two acquisitions.

In particular, the present invention achieves the aforementioned objective by exploiting, in addition to spectral and spatial information, also temporal information contained in the very frequent acquisitions made by the geostationary sensors, in order to detect even small fires (much smaller than the sensor spatial resolution), clouds and ground temperature with robustness and accuracy. In particular, the present invention is based on the combined use of a physical model of the radiative transfer process and a purely mathematical adaptive predictive algorithm to invert (solve) the RTM described by equations (1) and (5) exploiting a high number of acquisitions, much more than two. With respect to the technique based on equation (6), this inversion technique based on many acquisitions and on spectral and spatial information allows the physical parameters of interest to be estimated with higher accuracy and robustness.

BRIEF DESCRIPTION OF THE DRAWINGS

For a better understanding of the present invention, preferred embodiments, which are intended purely by way of example and are not to be construed as limiting, will now be described with reference to the attached drawings (all not to scale), wherein:

FIG. 1( a) shows schematically different contributions to the radiance acquired by a satellite sensor;

FIG. 1( b) shows the geometry of the down-welling thermal radiance emitted by the atmosphere;

FIG. 2( a) shows radiances of a region of the Earth with successive fire activities acquired by MSG/SEVIRI sensor during a day;

FIGS. 2( b) and 2(c) show estimated background temperature and fire extension, respectively, in the region of the Earth of FIG. 2( a), according to a first preferred embodiment of the invention;

FIGS. 3( a) and 3(b) show plots of radiances of a region of the Earth measured and predicted according to a second preferred embodiment of the invention;

FIG. 4 shows a table listing the fire detection probability corresponding to the overall burned area according to the present invention; and

FIG. 5 shows a table listing the false detection rate corresponding to the estimated fire extension according to the present invention.

DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS OF THE INVENTION

The following description is aimed at enabling a person skilled in the art to make and use the invention. Various modifications to the embodiments will be readily apparent to those skilled in the art, without departing from the scope of the present invention as claimed. Thus, the present invention is not intended to be limited to the embodiments shown, but is to be accorded the widest scope consistent with the principles and features disclosed herein and defined in the appended claims.

In order to improve the reliability of the estimation of the fire parameters, in the present invention a differential or, more precisely, a finite difference form of the Dozier RTM is conveniently used, where differences of the radiances between two acquisitions at close times are considered:

ΔR _(λ,t) ≡R _(λ,t) −R _(λ,t-Δt)=ε_(F,λ)·τ_(λ,t) ·B _(λ)(T _(F))·[f _(t) −f _(t-Δt)]++ε_(B,λ)·τ_(λ,t) ·[B _(λ)(T _(B,t))·(1−f _(t))−B _(λ)(T _(B,t-Δt))·(1−f _(t-Δt))]+R _(S,λ,t) −R _(S,λ,t-Δt)  (7)

In equation (7), the atmospheric term R_(A,λ) has been neglected because assumed unchanged within few acquisitions (order of tens of minutes). Moreover, since different bands are considered, the spectrally uncorrelated part of the disturbances is filtered out in the solution of (7). This equation can be considered for SEVIRI channels in the atmospheric windows in order to form a solvable equation system.

In order to remove the rough approximations affecting equation (6) and the assumptions on which it is based, several parameters have been introduced. Indeed, the ground temperatures T_(B,t) and T_(B,t-Δt) different times are considered different, the solar terms R_(S,λ,t) and R_(S,λ,t-Δt) are not neglected, and the fire and ground emissivities ε_(F,λ) and ε_(B,λ) are considered different.

Equation (7) has more unknown quantities than equation (6) and can be solved by means of a dynamic system approach. Indeed, the differential (or finite difference) Dozier RTM equation (7) can be considered as a dynamic system with the state variables f_(t) and T_(B,t):

$\begin{matrix} {{\begin{pmatrix} f_{t} \\ T_{B,t} \end{pmatrix} = {F\left( {\begin{pmatrix} f_{t - {\Delta \; t}} \\ T_{B,{t - {\Delta \; t}}} \end{pmatrix},T_{F},{\Delta \; R_{\lambda,t}},{\Delta \; R_{S,\lambda,t}},ɛ_{F,\lambda},ɛ_{B,\lambda},{\tau_{\lambda,t};\Lambda}} \right)}},} & (8) \end{matrix}$

where t and t−Δt denote two close acquisition times, ΔR_(λ,t)=R_(λ,t)−R_(λ,t-Δt), ΔR_(S,λ,t)=R_(S,λ,t)−R_(S,λ,t-Δt), and Λ is a set of wavelengths corresponding to the channels available in the atmospheric windows.

The dynamic equation (8) describes the time dynamic behavior of the Dozier RTM, where the quantities at the time t (f_(t) and T_(B,t)) can be estimated by means of the knowledge of the quantities calculated at the previous time t−Δt. In particular, for instance, the solar term can be calculated by means of equation (3) and emissivities as described in the aforementioned Early fire detection system based on multi-temporal images of geostationary and polar satellites, but other method/models can be considered. The transmittance, for instance, can be estimated by means of a parameterization similar to that in equation (2). Sensitivity analysis of equation (8) and test on real data has shown that the fire temperature T_(F) can be considered fixed (for instance at about 700° K) in order to increase the accuracy of the pixel fraction estimation. However, the fire temperature T_(F) could be considered as an unknown quantity, so further generalizing the dynamic equation (8).

The dynamic system approach expressed in equation (8) represents a physical model of the radiative transfer process.

The introduction of further unknown quantities, as shown in equation (7), and the use of the dynamic system approach shown in equation (8) allows the pixel fraction and background temperature to be estimated with higher accuracy and robustness with respect to the approach based on equation (6).

At each acquisition, at least two bands are needed to solve equation (8) and to estimate the state variables f_(t) and T_(B,t). If further channels are available, their exploitation makes the solution more accurate. The MIR channel (3.9 μm wavelength) is particularly sensitive to the presence of fires and it is used in every algorithm for fire detection from remote sensors. Unfortunately, in the SEVIRI sensor the MIR channel has a low saturation level that does not allow its use for monitoring large fires. Equation (8) can be solved even if the MIR channel cannot be used because saturated. Moreover, the solution of equation (8) is particularly robust to data gaps, occurring when some acquisitions are missing or cannot be used due to cloud coverage. In fact, the fire parameters can be still estimated by using sufficiently close acquisitions instead of two consecutive ones.

The accuracy of estimation of the pixel fraction f_(t) can be notably improved if oscillations with daily period are suppressed with a high-pass filter. This filtering is performed by removing from the last pixel fraction estimated by equation (8) the average of the previous consecutive N available pixel fractions (for instance 5 pixel fractions are enough):

$\begin{matrix} {{\hat{f}}_{t} = {f_{t} - {\frac{1}{N}{\sum\limits_{i = 1}^{N}f_{t - {\Delta \; t_{i}}}}}}} & (9) \end{matrix}$

The filtered pixel fraction {circumflex over (f)}_(t) is usefully exploited for fire detection. In fact, the accuracy of the filtered pixel fraction estimation obtained by using four bands of the SEVIRI/MSG sensor (1 MIR and 3 TIR channels) in the atmospheric windows is about 10⁻⁵.

FIG. 2( a) shows the SEVIRI radiances (brightness temperature of the middle infrared and thermal infrared bands 3.9 μm, 8.7 μm, 10.8 μm, 12 μm) acquired during about a day in a pixel with successive fire activities, while FIGS. 2( b) and 2(c) show the corresponding estimated background temperature and the estimated pixel fraction, respectively.

A fire is detected when the filtered pixel fraction {circumflex over (f)}_(t) is above a given threshold. A good compromise between probabilities of detection and of false alarm has been obtained with a threshold that corresponds to an active fire size of about 2000 m² using MSG/SEVIRI data.

As previously said, accuracy of the estimation of the fire parameters and detection and false alarm performances are significantly reduced when clouds, cloud borders, thin clouds, veils, fog (or the like) and low transmittance atmospheric conditions arise. A reliable cloud masking procedure is necessary to identify the acquisitions that are compatible with the physical model assumptions.

In order to recognize these unfavourable atmospheric conditions, an adaptive predictive algorithm is proposed which exploits the temporal information and correlations contained in a high number of acquisitions in order to detect clouds and fires.

The idea is that radiances vary slowly during a day and have an intrinsic periodicity of about one day, while the presence of clouds and/or fires adds high frequency oscillations, which also include data gaps. Quasi-periodicity of observed radiances during a clear-sky day without fires can be modeled by using few harmonic functions, and the model needs to be adapted in order to follow seasonal and other slow variations.

In particular, the adaptive predictive algorithm allows to realize an adaptive predictive model, wherein, exploiting SEVIRI/MSG data, a 24-hour spectral radiance sequence made up of a plurality of acquisitions at 96 acquisition times (every 15 minutes) is considered. The model is based on the estimation and filtering of the spectrum of the 24-hours spectral radiance sequence. The spectrum is evaluated by means of a Discrete Fourier Transform; then a low pass filter selects the harmonics of interest and, finally, the low-pass filtered signal is back-transformed. The adaptive predictive model is the back-transformed, low-pass filtered signal which is able to provide radiance predicted values.

In the following the adaptive predictive algorithm will be described step by step.

The first step consists of forming, for each pixel, a vector h_(λ) (h_(λ)(n), n=0, . . . , 95) containing a sequence of 96 time-consecutive spectral radiances R_(λ) relating to 24 hours and not necessarily acquired during the very same day. Should some acquisitions be unavailable, the relative spectral radiances may be computed by interpolation based on adjacent acquisitions.

In fact, for each pixel, several radiances acquired by the satellite sensor at different times and not affected by clouds or fires are needed.

In particular, the vector h_(λ) may be conveniently formed by considering several neighbouring days. A clouds-free acquisition for each vector element can be recognized as the one assuming the maximum spectral radiance value in a TIR band among those corresponding to the same vector element in the considered days. An analogous criterion is adopted to recognize fire-free acquisitions in a MIR band. A fire-free acquisition for each clouds-free vector element already selected can be recognized as the one assuming the minimum spectral radiance value in a MIR band among those corresponding to the same vector element in the considered days.

The second step of the adaptive predictive algorithm consists of computing the Discrete Fourier Transform (DFT) of the vector h_(λ), thus obtaining a vector H_(λ) with 96 elements defined as:

$\begin{matrix} {{{H_{\lambda}(k)} = {\sum\limits_{n = 0}^{N - 1}{{h_{\lambda}(n)}^{{- j}\frac{2\pi}{N}{nk}}}}},{k = 0},\ldots \mspace{14mu},95} & (10) \end{matrix}$

The third step consists of computing the adaptive predictive model associated to the λ μm band. The adaptive predictive model is a vector m_(λ) made up of 96 elements and defined as:

$\begin{matrix} {{{m_{\lambda}(n)} = {{\frac{1}{N}{H_{\lambda}(k)}} + {\frac{1}{N}{\sum\limits_{k = 1}^{A}\begin{bmatrix} {{{H_{\lambda}(k)}^{j\frac{2\pi}{N}{nk}}} +} \\ {{H_{\lambda}\left( {N - k} \right)}^{j\frac{2\pi}{N}{n{({N - k})}}}} \end{bmatrix}}}}},} & (11) \end{matrix}$

In particular, m_(λ) is the Inverse Discrete Fourier Transform of the low-pass filtered Discrete Fourier Transform of the vector h_(λ). The low-pass filter selects the continuous component (harmonic 0) and the first A harmonics (harmonic 1, 2, . . . , A) of the Discrete Fourier Transform of the vector h_(λ). For instance, with A=2, the model is reliable enough to fit a non-cloudy radiance trend during 24 hours.

Each n^(th) sample of the model m_(λ) contains the predicted value for the spectral radiance acquired at a given time of the day. For this reason, it is necessary to store the index n which points the element of m_(λ) relative to the last acquisition. In particular, assuming that the current acquisition at the time t is relative to the index n, the predicted radiance at the time t is given by the value m_(λ)( n):

$\begin{matrix} {R_{{PRD},\lambda} = {{m_{\lambda}\left( \overset{\_}{n} \right)} = {\frac{H_{\lambda}(0)}{N} + {\frac{1}{N}{\sum\limits_{k = 1}^{A}\begin{bmatrix} {{{H_{\lambda}(k)}^{j\frac{2\pi}{N}\overset{\_}{n}k}} +} \\ {{H_{\lambda}\left( {N - k} \right)}^{j\frac{2\pi}{N}{\overset{\_}{n}{({N - k})}}}} \end{bmatrix}}}}}} & (12) \end{matrix}$

The fourth step of the adaptive predictive algorithm is the cloud/fire detection which consists of checking if the predicted spectral radiance R_(PRD,λ) differs from the corresponding acquired spectral radiance R_(λ) using some thresholds. The cloud detection considers data acquired at 10.8 μm or at the 12.0 μm, while the fire detection considers data acquired at 3.9 μm.

In particular, the thresholds used in such detection procedure are three. The first one th_(DET,10.8) and the second one th_(DET,12) are used for the cloud detection allowing the comparison between the predicted values at 10.8 μm or at 12 μm and the spectral radiances acquired in such bands. The third threshold th_(DET,3.9) is used for the fire detection in a similar way in the band 3.9 μm.

More in detail, comparing the predicted value R_(PRD,λ) with the acquired spectral radiance R_(λ), a cloud is detected

-   -   if R_(PRD,10.8)−R_(10.8)>th_(DET,10.8); or     -   if R_(PRD,12)−R₁₂>th_(DET,12);         and a fire is detected, also if a cloud has been detected,     -   if R_(3.9)−R_(PRD,3.9)>th_(DET,3.9).

The adaptive predictive algorithm is able to correctly detect two different types of clouds: thick and thin clouds.

FIGS. 3( a) and 3(b) show plots of radiances measured and predicted according to the adaptive predictive algorithm during cloudy days. In particular, in FIG. 3( a) is shown the detection of clouds in absence of fire, while in FIG. 3( b) is shown the detection of both clouds and fires.

The fifth step of the adaptive predictive algorithm consists of updating the vector h_(λ) in order to follow seasonal and other slow variations of the radiances. If an acquisition is valid and it is not identified as a fire or a cloud, it is used to update the vector h_(λ). Also in this procedure three thresholds, th_(UPD,10.8), th_(UPD,12) and th_(UPD,3.9), are used to identify valid data to be used for updating the vector h_(λ).

In particular, acquired radiances R_(10.8), R₁₂ and R_(3.9) are considered valid if no anomalies are detected:

-   -   an anomaly at 10.8 μm is detected if         R_(PRD,10.8)−R_(10.8)>th_(UPD,10.8);     -   an anomaly at 12 μm is detected if R_(PRD,12)−R₁₂>th_(UPD,12);         and     -   an anomaly at 3.9 μm is detected if         R_(3.9)<R_(PRD,3.9)−th_(UPD,3.9) or if         R_(3.9)>R_(PRD,3.9)+th_(UPD,3.9).

If the current radiance R_(λ) is valid, it replaces in the vector h_(λ) the radiance contained in the position n (index relative to the current acquisition R_(λ)); otherwise the predicted radiance R_(PRD,λ) replaces in the vector h_(λ) the radiance value contained in the position n.

At each new acquisition, the pointer n is increased modulo 96 and the above described procedure is repeated starting from the second step of the adaptive predictive algorithm.

Note that the algorithm can be generalized and the tests on the validity of the radiances or on the presence of fires or clouds can be applied also to acquisitions at times before the last one, in order to refine previously obtained results based also on the new data as they become available.

The adaptive predictive model has high tracking capabilities even with long time acquisitions series affected by cloud coverage or missing data.

The adaptive predictive algorithm is used as a pre-processing step to determine when physical model can be applied. If no clouds are detected by the adaptive predictive algorithm in the analyzed pixel, along with the adaptive predictive model also the physical model can be exploited in order to detect fires and estimate the fire size (pixel fraction f_(t)) and the background temperature T_(B); otherwise, when clouds cover the analyzed pixel, only the adaptive predictive algorithm can still detect powerful fires. An example of this second case is shown in FIG. 3( b), where a fire is detected under thin clouds. The analysis of the 10.8 μm band highlights the presence of clouds that are correctly detected by the adaptive predictive algorithm. The oscillating results are due to the alternating occurrence of small clouds and clear sky. A big fire is detected by the adaptive predictive algorithm, though clouds are present and the physical model cannot be used.

A validation of the automatic fires detection method hereby disclosed has been performed using ground truth data relative to the Italian regions Sardinia, Lazio and Calabria, provided by the Italian Civil Protection Department. Two periods of 15 consecutive days in July 2004 and 15 consecutive days in August 2005 have been considered. These periods were characterized by not very cloudy days and a lot of occurrences of fires, in particular 286 fires occurred in areas greater than one hectare.

The validation has been carried out comparing the results obtained by the automatic fires detection method and the ground truth data. In particular, the fires have been classified in five groups depending on the size of the burned area as reported in the ground truth data. Concerning this, it is important to note that the automatic fires detection method detects the size of the flame frontline, which is the size of the real active fire and is much smaller than the overall burned area, while ground truth data reported only the overall burned area size.

For each burned area size class the probability of detection has been estimated as the detection rate: P_(d)=N_(TP)/(N_(TP)+N_(FN)), where N_(TP) is the number of true positive cases and N_(FN) the number of the false negative ones. The results are reported in the table shown in FIG. 4. As shown in FIG. 4, the fire detection probability increases with the increase of the overall burned area. The probability of false alarm has been estimated by means of the false alarm rate: P_(fa)=N_(FP)/(N_(FP)+N_(TN)), where N_(FP) is the number of false positive occurrences and N_(TN) is the number of the true negative ones.

The probability of false alarm does not depend on the fire size, in fact the estimation found in the performed validation has been P_(fa)=5.36·10⁻⁶, that means only 5 pixels per million erroneously detected as fires.

Another interesting index of the reliability of the automatic fires detection method is the false detection rate, defined as P_(fd)=N_(FP)/(N_(FP)+N_(TP)). This quantity depends on the estimated fire size, i.e. on the estimated pixel fraction (f_(t)). In FIG. 5 is shown a table reporting the false detection rate for different estimated fire sizes. As shown in FIG. 5, the false detection rate decreases with the increase of the estimated fire size.

From the foregoing, it may be immediately appreciated that the automatic fires detection method hereby disclosed allows to perform a reliable and robust fire detection both in presence and in absence of clouds. In particular, in absence of clouds covering the analyzed pixel, the method is able to reliably detect a fire and estimate the fire size and the background temperature, while, in presence of thin clouds, it is still able to detect fires. Moreover, the adaptive predictive algorithm is able to correctly detect the presence of thin and thick clouds in the analyzed pixel.

Finally, it is clear that numerous modifications and variants can be made to the present invention, all falling within the scope of the invention, as defined in the appended claims. 

1. A method for automatically detecting fires on Earth's surface using a satellite system, the method comprising: acquiring multi-spectral images of the Earth at different times, using a multi-spectral satellite sensor, each multi-spectral image being a collection of single-spectral images each associated with a respective wavelength (λ), and each single-spectral image being made up of pixels each indicative of a spectral radiance (R_(λ)) from a respective area of the Earth; providing a model relating the spectral radiances (R_(λ)) of each pixel in multi-spectral images acquired at different times and physical quantities representing thermodynamic phenomena occurring on the Earth's surface, in the Earth's atmosphere, and related to the Earth and the Sun relative positions; and computing, for at least one considered pixel, at a given time, at least the physical quantity in the model representing a possible fire on the Earth's surface, based on at least the spectral radiance (R_(λ,t)) at the given time, a previously acquired spectral radiance (R_(λ,t-Δt)), and the same physical quantity previously computed based on the same model.
 2. The method according to claim 1, wherein providing a model comprises: providing a dynamic multi-spectral equation system comprising at least one single-spectral equation associated with a respective wavelength (λ), and relating the spectral radiances (R_(λ)) of the considered pixels in single-spectral images acquired at different times to the physical quantities at the same times; and wherein computing the physical quantity in the model representing a possible fire on the Earth's surface comprises: computing at least the physical quantity representing a possible fire on the Earth's surface, based on the dynamic multi-spectral equation system.
 3. The method according to claim 2, wherein the single-spectral equations in the dynamic multi-spectral equation system are based on differences between equations relating the spectral radiances (R_(λ)) acquired at different times to the physical quantities at the same times.
 4. The method according to claim 2, wherein the computed physical quantities comprise a pixel fraction (f) corresponding to a fire extension and a background temperature (T_(B)).
 5. The method according to claim 4, wherein the physical quantities in the model further comprise one or more of the following quantities: a fire temperature (T_(F)); a Sun radiance (R_(S,λ)); a fire emissivity (ε_(F,λ)); a background emissivity (ε_(B,λ)); and an Earth's atmosphere transmittance (τ_(λ)).
 6. The method according to claim 4, further comprising: detecting a fire in an area of the Earth's surface whether the pixel fraction (f) computed for the corresponding pixel satisfies a given criterion.
 7. The method according to claim 6, wherein detecting a fire in an area of the Earth's surface whether the pixel fraction (f) computed for the corresponding pixel satisfies a given criterion includes: filtering the pixel fraction (f) computed for a pixel by subtracting from the computed pixel fraction (f) a value correlated to a low-pass filtered sequence of previously computed pixel fractions (f); and detecting a fire in an area of the Earth's surface whether the filtered pixel fraction ({circumflex over (f)}) computed for the corresponding pixel is higher than a given threshold.
 8. The method according to claim 5, wherein the dynamic multi-spectral equation system is: $\begin{pmatrix} f_{t} \\ T_{B,t} \end{pmatrix} = {F\left( {\begin{pmatrix} f_{t - {\Delta \; t}} \\ T_{B,{t - {\Delta \; t}}} \end{pmatrix},T_{F},{\Delta \; R_{\lambda,t}},{\Delta \; R_{S,\lambda,t}},ɛ_{F,\lambda},ɛ_{B,\lambda},{\tau_{\lambda,t};\Lambda}} \right)}$ wherein t and t−Δt denote two different times, λ denotes the wavelength, f_(t) and f_(t-Δt) denote the pixel fraction, T_(B,t) and T_(B,t-Δt) denote the background temperature, T_(F) denotes the fire temperature, ΔR_(λ,t)=R_(λ,t)−R_(λ,t−Δt) denotes a spectral radiance difference, ΔR_(S,λ,t)=R_(S,λ,t)−R_(S,λ,t-Δt) denotes a Sun radiance difference, ε_(F,λ) denotes the fire emissivity, ε_(B,λ) denotes the background emissivity, τ_(λ,t) denotes the Earth's atmosphere transmittance and Λ denotes the set of wavelengths (λ) considered; and wherein: the pixel fraction (f_(t)) and the background temperature (T_(B,t)) at the time t are obtained by solving the dynamic multi-spectral equation system; the spectral radiance difference (ΔR_(λ,t)) is directly computed from the acquired multi-spectral images; the fire temperature (T_(F)), the Sun radiance difference (ΔR_(S,λ,t)), the fire emissivity (ε_(F,λ)), the background emissivity (ε_(B,λ)) and the Earth's atmosphere transmittance (τ_(λ,t)) are computed independently; the pixel fraction (f_(t-Δt)) and the background temperature (T_(B,t-Δt)) at the time t−Δt are computed by solving the dynamic multi-spectral equation system at a previous iteration wherein at the first iteration the fire extension (f_(t-Δt)) and the background temperature (T_(B,t-Δt)) are set to respective values independently.
 9. The method according to claim 1, wherein providing a model comprises: computing an adaptive predictive model based on spectral radiances (R_(λ)) of pixels in multi-spectral images acquired at different times, the adaptive predictive model being able to predict spectral radiances (R_(PRD,λ)) of the considered pixels at a given time based on acquired spectral radiances (R_(λ)) and on previously computed predicted spectral radiances (R_(PRD,λ)); comparing the spectral radiance (R_(λ)) of the considered pixels at a given time with spectral radiances (R_(PRD,λ)) predicted for the same time by the adaptive predictive model; and detecting a fire in areas of the Earth's surface or an atmospheric phenomenon in areas of the Earth's atmosphere corresponding to the considered pixels based on an outcome of the comparison.
 10. The method according to claim 9, wherein computing an adaptive predictive model includes: for each considered pixel, computing a temporal harmonic analysis of the spectral radiances of the pixel related to different times, and acquired at, or predicted for the same wavelength (λ); and computing the adaptive predictive model based on computed low-frequency components of the spectral radiances so as to filter out high-frequency changes in the spectral radiances due to fires on Earth's surface or atmospheric phenomena.
 11. The method according to claim 10, wherein the temporal harmonic analysis is a Fourier analysis.
 12. The method according to claim 9, wherein detecting a fire in an area of the Earth's surface or an atmospheric phenomenon in an area of the Earth's atmosphere corresponding to a considered pixel comprises: detecting an atmospheric phenomenon if the acquired spectral radiance (R_(λ)) of the pixel is lower than the corresponding predicted spectral radiance (R_(PRD,λ)) of a given threshold value (th_(DET,λ)); and detecting a fire if the acquired spectral radiance (R_(λ)) of the pixel is higher than the corresponding predicted spectral radiance (R_(PRD,λ)) of a given threshold value (th_(DET,λ)).
 13. The method according to claim 12, wherein an atmospheric phenomenon is detected in an area of the Earth's atmosphere corresponding to a pixel by using the adaptive predictive model, the acquired spectral radiance of the pixel (R_(10.8)) and the threshold (th_(DET10.8)) associated with the 10.8 μm wavelength or by using the adaptive predictive model, the acquired spectral radiance of the pixel (R₁₂) and the threshold (th_(DET,12)) associated with the 12 μm wavelength; and wherein a fire is detected in an area of the Earth's surface corresponding to a pixel by using the adaptive predictive model, the acquired spectral radiance of the pixel (R_(3.9)) and the threshold (th_(DET,3.9)) associated with the 3.9 μm wavelength.
 14. The method according to claim 10, wherein computing, for each considered pixel, a temporal harmonic analysis includes: forming a vector (h_(λ)) containing the spectral radiances related to different times, and acquired at or predicted for the same wavelength (λ); and computing the temporal harmonic analysis of the vector (h_(λ)).
 15. The method according to claim 14, wherein forming a vector (h_(λ)) comprises: forming an initial vector (h_(λ)) containing spectral radiances (R_(λ)) of the considered pixel acquired at the same wavelength (λ) and at different times and not affected by atmospheric phenomena.
 16. The method according to claim 15, wherein forming an initial vector (h_(λ)) comprises: considering several neighbouring days; recognizing a clouds-free acquisition for each vector element as the one assuming the maximum spectral radiance value in the 10.8 or 12 μm wavelength band among those corresponding to the same vector element in the considered days; and recognizing a fire-free acquisition for each clouds-free vector element already selected as the one assuming the minimum spectral radiance value in 3.9 μm wavelength band among those corresponding to the same vector element in the considered days.
 17. The method according to claim 14, further comprising: classifying as valid or invalid acquired spectral radiances (R_(λ)) according to a given criterion; and wherein forming a vector (h_(λ)) further comprises: updating the vector (h_(λ)) with the acquired spectral radiances (R_(λ)) if they are classified as valid or with the corresponding predicted spectral radiances (R_(PRD,λ)) if the acquired spectral radiances (R_(λ)) are classified as invalid.
 18. The method according to claim 17, wherein classifying as valid or invalid acquired spectral radiances (R_(λ)) according to a given criterion comprises: classifying as valid the acquired spectral radiances (R₁₂) associated with the 12 μm wavelength if they are higher than the corresponding predicted spectral radiances (R_(PRD,12)) of a given threshold value (th_(UPD,12)); classifying as valid the acquired spectral radiances (R_(10.8)) associated with the 10.8 μm wavelength if they are higher than the corresponding predicted spectral radiances (R_(PRD,10.8)) of a given threshold value (th_(UPD,10.8)); and classifying as valid the acquired spectral radiances (R_(3.9)) associated with the 3.9 μm wavelength if they are in a range defined by the corresponding predicted spectral radiances (R_(PRD,3.9)) decreased and increased of a given threshold value (th_(UPD,3.9)).
 19. The method according to claim 17, wherein updating the vector (h_(λ)) comprises: replacing in the vector (h_(λ)) the spectral radiances at given times with the corresponding acquired spectral radiances (R_(λ)) classified as valid or with the corresponding predicted spectral radiances (R_(PRD,λ)).
 20. The method according to claim 9, further comprising: when no atmospheric phenomenon is detected, implementing also a method, wherein providing the model comprises: providing a dynamic multi-spectral equation system comprising at least one single-spectral equation associated with a respective wavelength (λ), and relating the spectral radiances (R_(λ)) of the considered pixels in single-spectral images acquired at different times to the physical quantities at the same times; and wherein computing the physical quantity in the model representing a possible fire on the Earth's surface comprises: computing at least the physical quantity representing a possible fire on the Earth's surface, based on the dynamic multi-spectral equation system, to automatically detect fires on the Earth's surface.
 21. A processing system configured to implement the method according to claim
 1. 22. A software loadable in a processing system and configured to implement, when run, the method according to claim
 1. 